Final Project workbook

# Load in data
hires_2019_clean <- read_csv('../clean_data/hires_2019_clean.csv')

Rows: 124446 Columns: 19
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr   (4): start_station_name, start_station_description, end_station_name, end_station_description
dbl  (11): duration, start_station_id, start_station_latitude, start_station_longitude, end_station_id, end_station_l...
date  (2): start_date, end_date
time  (2): start_time, end_time

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

rain_2019_clean <- read_csv('../clean_data/rain_2019_clean.csv')

Rows: 365 Columns: 2
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl  (1): rainfall_mm
date (1): start_date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(hires_2019_clean)

Rows: 124,446
Columns: 19
$ start_date                <date> 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, …
$ start_time                <time> 05:58:47, 06:02:13, 06:14:09, 06:25:08, 06:54:31, 07:31:24, 08:47:46, 08:48:06, 09:0…
$ end_date                  <date> 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, 2019-01-01, …
$ end_time                  <time> 07:16:14, 07:21:17, 07:21:57, 06:36:22, 06:59:52, 07:42:31, 09:00:57, 09:36:29, 09:3…
$ duration                  <dbl> 77.43, 79.07, 67.80, 11.23, 5.33, 11.12, 13.18, 48.37, 24.80, 81.05, 47.33, 15.42, 58…
$ start_station_id          <dbl> 259, 259, 251, 358, 251, 261, 259, 257, 358, 284, 257, 257, 257, 257, 257, 253, 246, …
$ start_station_name        <chr> "St Andrew Square", "St Andrew Square", "Waverley Station", "Leith Walk", "Waverley S…
$ start_station_description <chr> "North East corner", "North East corner", "Platform level near Calton Road Exit", "Le…
$ start_station_latitude    <dbl> 55.95475, 55.95475, 55.95264, 55.96504, 55.95264, 55.96093, 55.95475, 55.95700, 55.96…
$ start_station_longitude   <dbl> -3.192774, -3.192774, -3.187527, -3.176686, -3.187527, -3.181005, -3.192774, -3.11688…
$ end_station_id            <dbl> 285, 285, 285, 284, 261, 265, 246, 344, 257, 257, 284, 257, 225, 256, 256, 246, 247, …
$ end_station_name          <chr> "Ocean Terminal", "Ocean Terminal", "Ocean Terminal", "Leith Links", "Brunswick Place…
$ end_station_description   <chr> "Debenhams / Puregym entrance", "Debenhams / Puregym entrance", "Debenhams / Puregym …
$ end_station_latitude      <dbl> 55.98129, 55.98129, 55.98129, 55.96920, 55.96093, 55.93981, 55.93900, 55.93640, 55.95…
$ end_station_longitude     <dbl> -3.176351, -3.176351, -3.176351, -3.166885, -3.181005, -3.182739, -3.173924, -3.19425…
$ rainfall_mm               <dbl> 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0…
$ start_elevation           <dbl> 76.04, 76.04, 57.10, 20.96, 57.10, 33.11, 76.04, 4.30, 20.96, 15.45, 4.30, 4.30, 4.30…
$ end_elevation             <dbl> 7.30, 7.30, 7.30, 15.45, 33.11, 80.45, 88.97, 94.55, 4.30, 4.30, 15.45, 4.30, 53.83, …
$ elevation_diff            <dbl> -68.74, -68.74, -49.80, -5.51, -23.99, 47.34, 12.93, 90.25, -16.66, -11.15, 11.15, 0.…

# visualising use by month
hires_2019_clean %>% 
  mutate(month = month(start_date, label = TRUE), .before = 1) %>% 
  group_by(month) %>% 
  summarise(count = n()) %>%
  ggplot() +
  aes(x = month, y = count) +
  geom_col(fill = "#F27F1B") +
  geom_text(aes(label = count), vjust = 2, colour = "white", size = 3.5) +
  labs(x = "\nMonth",
       y = "Number of journeys\n",
       title = "Total number of bike journeys by month (2019)",
       subtitle = "Total = 124446 journeys\n\n") +
  scale_y_discrete(expand = c(0,1), limits = c(0, 5000, 10000, 15000)) +
  theme_minimal() +
  theme(title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

Mapping out stations

# create smaller dataset for map
map_data <- hires_2019_clean %>% 
  select(start_station_longitude, start_station_latitude, start_station_id, start_station_name, start_elevation) %>% 
  distinct(start_station_id, .keep_all = TRUE)

# create custom icon for bike hire stations
bike  <-  makeIcon("../www/bicycle-outline.png")

# map station data using leaflet app
map_data %>% 
  leaflet() %>% 
  addProviderTiles(providers$CartoDB.Positron) %>% 
  addTiles(providers$CartoDB.Positron) %>%
  addMarkers(~start_station_longitude, ~start_station_latitude, 
             icon = ~bike,
             #clusterOptions = markerClusterOptions(),
             popup = ~paste0("Station ID: ", start_station_id,
                              "<br>Name: ", start_station_name,
                              "<br>Elevation (m): ", start_elevation))

  # addCircles(lng = ~start_station_longitude,
  #            lat = ~start_station_latitude,
  #            color = "#F27F1B",
  #            popup = ~paste0("Station ID: ", start_station_id,
  #                            "<br>Name: ", start_station_name,
  #                            "<br>Elevation (m): ", start_elevation))

Most popular stations

# top 10 start and end stations
hires_2019_clean %>% 
  select(start_station_id, start_station_name, start_station_description, start_elevation) %>% 
  count(start_station_id, start_station_name, start_station_description, start_elevation) %>% 
  arrange(desc(n)) %>% 
  head(10)

hires_2019_clean %>% 
  select(end_station_id, end_station_name, end_station_description, end_elevation) %>% 
  count(end_station_id, end_station_name, end_station_description, end_elevation) %>% 
  arrange(desc(n)) %>%
  head(10)

hires_2019_clean %>% 
  select(start_station_id, end_station_id) %>% 
  filter(start_station_id == c("248", "259", "265", "289", "257", "262", "249", "247", "171", "183"), 
         end_station_id == c("262", "257", "250", "265", "248", "358", "259", "183", "171", "258")) %>% 
  chordDiagram(scale = TRUE)

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘select’ for signature ‘"spec_tbl_df"’

hires_2019_clean %>% 
  select(start_station_id, end_station_id) %>%
  filter(start_station_id == c("248", "259", "265", "289", "257", "262", "249", "247", "171", "183"), 
         end_station_id == c("262", "257", "250", "265", "248", "358", "259", "183", "171", "258")) %>% 
  ggraph(layout = "linear") + 
  geom_edge_arc(edge_colour = "black", edge_alpha = 0.3, edge_width = 0.2) +
  geom_node_point(color = "#F27F1B", size = 10) +
  geom_node_text(aes(label = name), repel = FALSE, size = 3, nudge_y = 0, colour = "white") +
  labs(title = "Journeys between most popular stations\n\n") +
  theme_void() +
  theme(title = element_text(size = 12),
        legend.position = "none",
        plot.margin = unit(rep(1, 4), "cm"))

Journeys based on elevation

68622 downhill journeys 43774 uphill journeys 12050 flat journeys

12050 / 124446 * 100

[1] 9.682915

Journeys based on weather

According to www.statista.com “A rainday is when one millimetre or more of rain occurs in a day.” Based on this I have classified rainy days as days where rainfall as recorded by www.power.larc.nasa.gov/ was equal to or over 1mm in Edinburgh.

Hypothesis Test

# calculate the p-value and compare to α value
p_value <- null_distribution %>%
  get_p_value(obs_stat = observed_stat, direction = "left")

Warning: Please be cautious in reporting a p-value of 0. This result is an approximation based on the number of `reps` chosen in the `generate()` step. See `?get_p_value()` for more information.

# hypothesis test
# α = 0.05
# H0 - Rain makes no difference to bike hires
# H1 - Rain has an impact on bike hires
# Type of test is one-sample proportion, right-sided.
# 
# p value = 0 
# 
# The p-value being near to 0 means it is less than α, and so we should reject H0 based on our data. The result is statistically significant in favour of H1, that rain may impact bike hire numbers.

# rain_day <- hires_2019_clean %>%
#   group_by(start_date) %>% 
#   mutate(rain_day = rainfall_mm >= 1) %>% 
#   select(start_date, rain_day) %>% 
#   count(start_date, rain_day)
# 
# null_distribution <- rain_day %>%
#   specify(response = rain_day, success = "TRUE") %>%
#   hypothesize(null = "point", p = 0.05) %>%
#   generate(reps = 10000, type = "draw") %>%
#   calculate(stat = "prop")
# 
# obs_stat <- rain_day %>%
#   specify(response = rain_day, success = "TRUE") %>%
#   calculate(stat = "prop")
# 
# null_distribution %>%
#   visualise() +
#   shade_p_value(direction = "right", obs_stat = obs_stat)
# 
# null_distribution %>%
#   get_p_value(direction = "right", obs_stat = obs_stat)

# bootstrapped hypothesis
# α = 0.05
# H0 - Rain makes no difference to bike hires
# H1 - Rain has an impact on bike hires

# rain_day_flag <- rain_day %>%
#   mutate(rain_day_flag = if_else(rain_day == "TRUE", 1, 0))
# 
# null_distribution <- rain_day_flag %>%
#   specify(response = rain_day_flag) %>%
#   hypothesize(null = "point", mu = 0.05) %>%
#   generate(reps = 10000, type = "bootstrap") %>%
#   calculate(stat = "mean")
# 
# null_distribution %>%
#   visualise() +
#   shade_p_value(direction = "right", obs_stat = obs_stat)
# 
# null_distribution %>%
#   get_p_value(direction = "right", obs_stat = obs_stat)

Independence Hypothesis Test On Rainy day data

Does a rainy day have an effect on bike hires

α = 0.05 H0 - Rain makes no difference to bike hires H1 - Rain has an impact on bike hires

The p-value ≤ α, so I reject the null hypothesis H0 in favour of the alternative hypothesis H1, that rain may have an effect on bike hires.

Rainy Days in 2019 = 181 Rainy journeys = 58073 (46.7%) Non rainy journeys = 66373 (53.3%)

66373 / 124446 * 100

[1] 53.33478

Time of Day analysis

Gym Bunnies - 5am to 7am: 2419 17.54 mins
Morning Commuters - 7am to 9am: 12452 17.87 mins
Day Trippers - 9am to 5pm: 66880 29.75 mins
Homeward Bounders - 5pm to 6.30pm: 14961 23.86 mins
Evening Movers - 6.30pm to 10pm: 16815 23.50 mins
Pub Ponies - 10pm to 5am: 10919 23.07 mins

Average journey time overall is 26.18 minutes

Most popular hour is from 4pm to 5pm 12279 journeys average time 24.27 mins.

hires_2019_clean %>% 
  select(start_time) %>%
  group_by(start_time) %>% 
  count() %>%
  ggplot() +
  geom_histogram(aes(start_time), col = "white", fill = "#F27F1B") +
  labs(x = "\nTime of Day",
       y = "Number of Journeys\n",
       title = "Start time of journeys") +
  theme_minimal() +
  theme(title = element_text(size = 12),
        axis.text = element_text(size = 10),
        axis.title.x = element_text(size = 10),
        axis.title.y = element_text(size = 10),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "black"))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

hires_2019_clean %>% 
  mutate(start_hour = hour(start_time), .before = 2) %>% 
  group_by(start_hour) %>% 
  count(start_hour) %>% 
  ggplot() +
  #geom_rect(aes(xmin=0, xmax=7, ymin=0, ymax=Inf), fill = "#F27F1B", alpha = 0.1) +
  # annotate("rect", xmin = 0, xmax = 5, ymin = 0, ymax = Inf, fill = "#5A4FF5", alpha = .9) +
  # annotate("rect", xmin = 5, xmax = 7, ymin = 0, ymax = Inf, fill = "#FFFB85", alpha = .7) +
  # annotate("rect", xmin = 7, xmax = 9, ymin = 0, ymax = Inf, fill = "#FFE930", alpha = .6) +
  # annotate("rect", xmin = 9, xmax = 17, ymin = 0, ymax = Inf, fill = "#F27F1B", alpha = .4) +
  # annotate("rect", xmin = 17, xmax = 18.5, ymin = 0, ymax = Inf, fill = "#DE6DCC", alpha = .4) +
  # annotate("rect", xmin = 18.5, xmax = 22, ymin = 0, ymax = Inf, fill = "#9D84F5", alpha = .7) +
  # annotate("rect", xmin = 22, xmax = 23, ymin = 0, ymax = Inf, fill = "#5A4FF5", alpha = .9) +
  geom_point(aes(start_hour, n), color = "#F27F1B") +
  geom_line(aes(start_hour, n), color = "#F27F1B", alpha = 0.5) +
    labs(x = "\nStart hour",
       y = "Number of Journeys\n",
       title = "START TIME OF JOURNEYS\n") +
  scale_x_discrete(expand = c(0,.5), limits = c(0:23)) +
  theme_minimal() +
  theme(title = element_text(size = 12, family="Proxima Nova", face = "bold"),
        axis.text = element_text(size = 10, family="Proxima Nova", face = "bold", colour = "#8C8C8C"),
        axis.title.x = element_text(size = 10, colour = "#8C8C8C"),
        axis.title.y = element_text(size = 10, colour = "#8C8C8C"),
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(), 
        axis.line = element_line(colour = "#8C8C8C"))

Warning: Continuous limits supplied to discrete scale.
Did you mean `limits = factor(...)` or `scale_*_continuous()`?

# trying a raster image of elevation

edi <- st_read('../raw_data/Edinburgh_Ward_Boundaries/Edinburgh_Ward_Boundaries.shp', "Edinburgh_Ward_Boundaries")

st_union(edi)

ggplot() +
  geom_raster(data = new_elevation_raster, aes(x = x, y = y)) +
  geom_sf(data = edi, color = "red") +
  coord_sf() +
  scale_fill_viridis_c() +
  labs(title = " ", x = " ", y = " ", fill = " ")

Warning: Raster pixels are placed at uneven vertical intervals and will be shifted. Consider using geom_tile() instead.

plot(elevation_raster, 
     main = "Elevation levels across Edinburgh in metres",
     key = "map_data")

Warning in plot.window(...) : "key" is not a graphical parameter
Warning in plot.xy(xy, type, ...) : "key" is not a graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...) :
  "key" is not a graphical parameter
Warning in axis(side = side, at = at, labels = labels, ...) :
  "key" is not a graphical parameter
Warning in box(...) : "key" is not a graphical parameter
Warning in title(...) : "key" is not a graphical parameter